/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#include "qfunc.h"
#include "polytomousrasch.c"

#include <string.h>

double
qfunc (double                           kappa,
       const double*                    logalpha,
       unsigned int                     n_workers,
       unsigned int                     n_labels,
       const double*                    tau,
       unsigned int                     n_images,
       const double*                    gamma,
       const double*                    pz,
       const double*                    priorz,
       const Rating*                    ratings,
       unsigned int                     n_ratings,
       const Distribution*              priorkappa,
       const ParameterizedDistribution* priorlogalpha,
       const ParameterizedDistribution* priortau,
       Distribution*const*              priorgamma)
{
  double q = 0;

  /* TODO: This first loop computes a constant as far as the minimization
   * routine is concerned.  Maybe it can be eliminated?  (Maybe why not,
   * the termination condition is a function of the relative change in q).
   */

  for (unsigned int i = 0; i < n_images; ++i)
    {
      for (unsigned int j = 0; j < n_labels; ++j)
        {
          q += pz[j + i * n_labels] * log (priorz[j]);
        }
    }

  for (unsigned int i = 0; i < n_ratings; ++i)
    {
      double alpha = exp (logalpha[ratings[i].image]);

      for (unsigned int j = 0; j < n_labels; ++j)
        {
#define taufunc(_k)                                             \
  ({                                                            \
     alpha * tau[ratings[i].worker * (n_labels - 1) + _k];      \
  })
          q += pz[j + ratings[i].image * n_labels] * 
               log_polytomous_rasch (ratings[i].rating,
                                     alpha * j,
                                     taufunc,
                                     n_labels - 1);
        }
#undef taufunc
    }

  q += distribution_logpdf (priorkappa, kappa);

  for (unsigned int i = 0; i < n_images; ++i)
    {
      q += parameterized_distribution_logpdf (priorlogalpha,
                                              kappa,
                                              logalpha[i]);
    }

  for (unsigned int i = 0; i < n_workers; ++i)
    {
      for (unsigned int j = 0; j + 1 < n_labels; ++j)
        {
          q += parameterized_distribution_logpdf (priortau,
                                                  gamma[j],
                                                  tau[i * (n_labels - 1) + j]);
        }
    }

  for (unsigned int j = 0; j + 1 < n_labels; ++j)
    {
      q += distribution_logpdf (priorgamma[j], gamma[j]);
    }

  return q;
}

void
dqfunc (double*                         dkappa,
        double*                         dlogalpha,
        double*                         dtau,
        double*                         dgamma,
        double                          kappa,
        const double*                   logalpha,
        unsigned int                    n_workers,
        unsigned int                    n_labels,
        const double*                   tau,
        unsigned int                    n_images,
        const double*                   gamma,
        const double*                   pz,
        const Rating*                   ratings,
        unsigned int                    n_ratings,
        const Distribution*             priorkappa,
        const ParameterizedDistribution* priorlogalpha,
        const ParameterizedDistribution* priortau,
        Distribution* const*            priorgamma)
{
  *dkappa = 0;
  memset (dlogalpha, 0, n_images * sizeof (*dlogalpha));

  if (n_labels > 0)
    {
      memset (dtau, 0, n_workers * (n_labels - 1) * sizeof (*dtau));
      memset (dgamma, 0, (n_labels - 1) * sizeof (*dgamma));
    }

  for (unsigned int i = 0; i < n_ratings; ++i)
    {
      double alpha = exp (logalpha[ratings[i].image]);

      for (unsigned int j = 0; j < n_labels; ++j)
        {
          double tmpdtau[n_labels - 1];
#define taufunc(_k)                                             \
  ({                                                            \
     alpha * tau[ratings[i].worker * (n_labels - 1) + _k];      \
  })
          dlogalpha[ratings[i].image] += 
            pz[j + ratings[i].image * n_labels] * 
            alpha * 
            j * 
            dlog_polytomous_rasch_dbeta (ratings[i].rating,
                                         alpha * j,
                                         taufunc,
                                         n_labels - 1);

          dlog_polytomous_rasch_dtau (tmpdtau,
                                      ratings[i].rating,
                                      alpha * j,
                                      taufunc,
                                      n_labels - 1);
#undef taufunc

          for (unsigned int k = 0; k + 1 < n_labels; ++k)
            {
              dtau[ratings[i].worker * (n_labels - 1) + k] +=
                pz[j + ratings[i].image * n_labels] * 
                alpha * 
                tmpdtau[k];
              dlogalpha[ratings[i].image] += 
                pz[j + ratings[i].image * n_labels] * 
                alpha *
                tau[ratings[i].worker * (n_labels - 1) + k] * 
                tmpdtau[k];
            }
        }
    }

  *dkappa += distribution_dlogpdf (priorkappa, kappa);

  for (unsigned int i = 0; i < n_images; ++i)
    {
      *dkappa +=
        parameterized_distribution_d0logpdf (priorlogalpha,
                                             kappa,
                                             logalpha[i]);
      dlogalpha[i] += 
        parameterized_distribution_d1logpdf (priorlogalpha,
                                             kappa,
                                             logalpha[i]);
    }

  for (unsigned int i = 0; i < n_workers; ++i)
    {
      for (unsigned int j = 0; j + 1 < n_labels; ++j)
        {
          dgamma[j] += 
            parameterized_distribution_d0logpdf (priortau,
                                                 gamma[j],
                                                 tau[i * (n_labels - 1) + j]);
          dtau[i * (n_labels - 1) + j] +=
            parameterized_distribution_d1logpdf (priortau,
                                                 gamma[j],
                                                 tau[i * (n_labels - 1) + j]);
        }
    }

  for (unsigned int j = 0; j + 1 < n_labels; ++j)
    {
      dgamma[j] += distribution_dlogpdf (priorgamma[j], gamma[j]);
    }
}
